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INTRODUCTION 

In  the  theory  of  linear  differential  equations  of  the 
elliptic  type,  boundary-value  problems  play  a  basic  role.   Such 
problems  involve  the  determination  of  the  solution  of  a  given 
differential  equation  within  a  given  domain  on  the  boundary  of 
which  given  conditions  must  be  fulfilled. 

Partial  differential  equations  in  two  variables  involve  a 
domain  of  application  in  a  two-dimensional  plane.   Boundary 
values  may  be  assigned  in  various  ways  along  curves  in  the  plane 
to  determine  the  solution  elsewhere  in  the  plane.   For  second- 
order  equations,  a  classification  into  three  basic  types  is 
useful  in  distinguishing  boundary  value  problems. 

Consider  a  second  order  partial  differential  equation  of 
the  form 


[1]    a(x,y)<j>xx+  b(x,y)<t>   +  0<x,y)6  +  f  (x,y  .*x>*y  3  =  0. 

i.e.,  one  that  is  linear  in  the  second-order  partial  derivatives. 

It  is  possible  to  classify  this  equation  on  the  basis  of  the 

2 
expression  b   -  tac,  where  a,  b  and  c  are  the  coefficients  of 

the  second  order  terms  of  the  equation.   The  classification  is 

as  follows : 

2 

1.  Parabolic  Type:   (if  b   -  4ac  =  0  at  all  points  of  the 

region  considered).   An  example  of  this  type  is  the  one-dimen- 
sional heat-flow  equation. 

2 

2.  Hyperbolic  Type:   (if  b   -  4ac  >  0  for  all  points  of 

the  region).   A  basic  problem  of  this  type  is  the  one-dimensional 


wave  equation. 

3.   Elliptic  Type:   (if  b   -  4ac  <  0  at  all  points  of  the 
region).   The  most  fundamental  equation  of  this  type  is  Laplace's 
equation.   The  elliptic  equation  will  be  considered  in  this 
report.   If,  in  the  elliptic  type  equation,,  b  and  f  are  both 
zero  and  a  =  c  =  1,  the  equation  is  called  Laplace's  equation  in 
two  dimensions  and  its  solution  is  referred  to  as  a  harmonic 
function. 

Any  elliptic  equation  with  f  =  0  can  be  reduced  to  Laplace's 
equation.   First  of  all,  a  rotation  of  axes  can  be  used  to 
remove  the  mixed  partial  term.   A  common  transformation  to 
accomplish  this  is 

x  ■  Xcos  <j>  -  Ysin  <j> ,      y  =  Xsin  $  +  Ycos  <t> 

2 

where  <j>  is  the  angle  of  rotation. 

The  general  equation  then  reduces  to  the  form 

a'(X,Y)<Cxx  +  c'(X,YHYV  =  0. 

Then  a  change  of  scale  can  be  used  to  make  the  coefficients  of 
the  second  partials  equal. 

Boundary-value  problems  of  this  type  involve  specifying 
values  of  the  function  or  of  its  partial  derivatives  of  first 


1.  A  solution  of  Laplace's  equation  (in  any  number  of 
dimensions)  is  called  a  harmonic  function. 

2.  Derived  in  Calculus ,  Tom.  M.  Apostol,  Blaisdell 
Publishing  Company,  New  York,  London,  1961,  pp.  318-320. 


order  along  the  closed  boundary  of  the  region.   Then  solutions 
are  determined  for  all  interior  points  of  the  region.   A  typical 
boundary  value  problem  for  elliptic  equations  is  the  Dirichlet 
problem  which  may  be  defined  as  follows: 

DEFINITION:   (Dirichlet  Problem):   "Let  f  be  a  continuous  function 
prescribed  on  the  boundary  B  of  some  finite  region  G  of  the  space 
(x-,,  x,,  ...  ,  x  ).   We  are  to  find  a  function  $(x^,    x2 ,  ...  ,  xn) 
which  is  harmonic  in  the  interior  of  G  and  equal  to  f  on  B." 

In  general,  numerical  methods  of  solution  of  all  Dirichlet 
boundary-value  problems  depend  on  the  use  of  a  grid  of  points  in 
the  domain.   At  the  nodes  of  this  grid,  approximations  to  the 
solution  are  determined.   With  an  infinitely  fine  mesh,  the 
solution  would  converge  to  the  solution  of  the  partial  differen- 
tial.  Various  coordinate  meshes  may  be  used  in  the  two  dimen- 
sional system;  however,  the  configuration  in  this  report  is 
restricted  to  a  square  lattice; 

Laplace's  equation  in  two  dimensions, 

2      2 

9  <t>  +  3  <fr  _  n 

2      2  "   ' 
Sx     3y 

2  2 

is  sometimes  written  7  0=0.   The  expression  v  ^  >  whether 

referring  to  a  rectangular  coordinate  system  as  above ,  in  two 

or  more  dimensions,  or  the  corresponding  expression  in  some  other 

system  of  coordinates,  is  called  the  Laplacian. 

In  a  square  lattice  the  two-dimensional  Laplacian  can  be 

approximated  by  replacing  the  second  derivatives  by  their 


expressions  in  terms  of  finite  differences.   It  is  customary  to 
use  only  the  approximation  to  the  Laplacian  obtained  by  neglect- 
ing the  fourth  and  higher  differences.   If  a  single  subscript  is 
used  to  designate  the  relative  position  of  the  <j> '  s  in  the  lattice, 
the  expression  to  evaluate  the  potential  function  <t>  at  each 
interior  point  looks  simpler,  compared  with  using  double  sub- 
scripts.  However,  a  system  must  be  defined  for  a  consistent 

numbering  pattern  of  the  lattice  points.   Then  the  Laplacian, 

3 
expressed  in  terms  of  finite  differences,  may  be  written  as 

'^0  =  i/h2  (^  +  *2  +  <p3  +  <s>K   -  m  <e0), 

where  <t>,  ,  <)>-,  ♦«,  and  <K  are  the  mesh  points  h  units  from  41-. 
In  symbolic  notation 


'h*0  '  1/h'   <  *3      -»  ' 


This  report  is  based  on  the  Liebmann  method  of  solving 
Laplace's  equation.   In  general,  the  Liebmann  method  is  a 
process  for  evaluating  the  potential  function  at  an  interior 
point  in  terms  of  its  neighboring  points.   The  method  and  its 
extensions  are  illustrated  in  the  following  sections. 


3.   Kaiser  S.  Kunz,  Numerical  Analysis ,  McGraw-Hill  Book 
Company,  Inc.,  New  York,  London,  Toronto,  1957,  pp.  278-79. 


LIEBMANN  METHOD  -  TWO  DIMENSIONAL 

The  method  of  Liebmann  can  best  be  illustrated  by  a  simple 
example,  derived  with  the  use  of  Laplace's  equation.   Assume 
that  the  potential  * ,  satisfying  Laplace's  equation,  is  required 
for  the  region  contained  within  a  given  square  boundary  (see 
figure  1).   The  potential  <f>  is  zero  on  three  sides  of  the  square. 
On  the  fourth  side  <j>  .«'.  1  at  the  middle  and  drops  off  linearly 
to  zero  at  the  corners . 
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In  this  example,  a  square  net  of  25  lattice  points  is  used 
to  approximate  the  behavior  of  <(>.   The  potential  at  the  interior 
lattice  points  can  be  designated  by  <|>  where  s  =  1,  2,  ...  ,9; 
the  potential  at  the  boundary  points  by  <j>  ,  t  =  10,  11,  ...  ,  25. 
The  problem  then  is  to  determine  a  numerical  solution  of  the  <J>g 
in  terms  of  the  $. . 

It  has  already  been  shown  that  Laplace's  equation  can  be 
approximated  by  the  partial  difference  equation  indicated  sym- 
bolically by 

V2,),   =  l/h2  <  1     -4  0, 


hYs 


1 


where  h  is  the  net  spacing.   Consider  the  potential  at  interior 
point  $. .   The  difference  equation  which  expresses  the  potential 
for  this  point  is  written 

l/h2    (<|>1;L  +   4>25   +   <j>2   +   <t>4   -   4   <j>1)    =   0. 

We  can  solve  this  equation  for  4  <j>,  and  write  it  as 

4  *1  =  *ll  +  *25  *  *2  +  U    ' 

or  write  it  in  the  following  form  which  allows  only  the  boundary 
points  on  the  right  hand  side  of  the  equation: 

[3]  -4  #j  +  <t>2  +  'c^  =  -  (*X1  +  *2g) 

Equation  3  gives  the  potential  at  the  first  interior  point 
$   =  1.   As  we  consider  the  remaining  interior  points, 


s=2,3,  ...  ,9,  we  have  eight  additional  equations. 

*1  -14  *2  +  *3  +  *5  =  "  *12 

*2  "4  *3  +  *6  =  "  (*13  +  W 

♦i  'H  K  +  h  +  h  '-  ~   *24 

*2  +  \    "4  *5  +  *6  +  *8  =  ° 

*3  +  *5  "4  *6  +  *9  =  "  *16 

\    "4  *7  +  *8  =  "  (*21  +  *23) 

*5  +  *7  "4  *s'+  *9  =  "  *20 

*6  +  *8  ~H    *9  =  "  (*17  +  *19) 
Notice  that  the  right  hand  sides  of  these  equations  are 
known  as  they  represent  values  at  the  boundary  points  which  are 
constant.   We  have  a  system  of  nine  equations  in  as  many  unknowns 
which  can  be  solved  by  standard  methods.   However,  if  the  system 
should  be  expanded  and  a  large  number  of  interior  points  con- 
sidered, the  method  would  become  an  exercise  in  "busy  work". 
For  this  reason  we  look  to  other  methods. 

Consider  the  set  of  linear  equations :   equation  3  and  the 
eight  additional  equations  for  this  problem.   The  matrix  of  the 
coefficients  for  this  system  of  equations  has  -4  as  the  elements 
of  its  principal  diagonal.   Notice  that  all  of  the  other  elements 
are  either  0  or  1.   Hence,  the  diagonal  elements  are  dominant. 
In  fact,  the  absolute  value  of  any  element  on  the  diagonal  is  at 
least  as  large  as  the  sum  of  the  remaining  elements  in  that  row 
or  column.   The  nature  of  this  coefficient  matrix  indicates  that 


we  can  safely  look  to  some  method  of  successive  approximations 
for  our  solution.    One  such  method  is  introduced  by  Liebmann. 

If  the  symbolic  difference  equation  expressed  above  is 
rewritten,  we  have ,■  symbolically , 


[4]  <j>a  =  1/4  \   : 


This  equation  states  that  $      is  just  an  average  of  the  values  of 
cj>  at  the  four  neighboring  points.   Neighboring  points  are  defined 
as  those  points  a  distance  h  away  from  <j> .   Liebmann 's  method 
consists  of  improving  the  value  initially  guessed  for  $  by 
repeated  application  of  this  process  over  the  set  of  points. 

One  passes  from  point  to  point  in  the  lattice  replacing 
the  previous  values  of  <f  at  each  point  by  the  average  of  the  cj> '  s 
for  the  four  closest  neighboring  points.   All  points  in  the 
interior  of  the  region  must  be  included  in  the  process  ,  although 
it  does  not  matter  what  sort  of  a  pattern  is  originally  set  up 
for  each  region. 

As  illustrated  in  figure  1,  one  possible  pattern  for 
numbering  the  lattice  points  is  to  begin  at  the  upper  left 
interior  point.   Number  consecutively  along  the  rows,  not  includ- 
ing any  point  falling  on  the  boundary  of  the  lattice.   The 


4.  Ibid. ,  pg.  296. 

5.  H.  Liebmann,  Die  Angenaherte  Ermittlung  Harmonischer 
Funktionen  und  Konformer  Abbildungen,  Sitzber,  Math.  Physik,  kl. 
Boyer.  Alcad .   Wiss.  Munchen,  47:   385-416  (1918). 


potential  is  evaluated  for  $. .   Then  the  potential  is  evaluated 
for  the  interior,  point  of  next  higher  index.   The  process  con- 
tinues in  the  defined  order  throughout  the  entire  interior  of  the 
lattice.   When  the  potential  at  all  points  has  been  improved  once, 
these  new  values  are  substituted  for  the  values  of  the  previous 
iteration  and  the  process  begins  again  at  <j>,  . 

The  process  converges  since  each  application  of  the  method 
averages  the  errors  of  the  four  neighboring  points .   If  the 
errors  are  not  all  of  the  same  sign,  they  will  tend  to  average 
out  and  <f  will  be  improved.   In  several  cases,  one  or  more  of 
the  four  neighboring  points  is  actually  a  boundary  point  which 
has  no  error.   Although  the  process  converges,  in  some  problems 
the  convergence  is  rather  slow. 

The  Liebmann  process  is  well  suited  to  a  large  scale 
computer:   however,  the  requirement  that  values  for  all  <j> '  s  in 
the  lattice  must  be  made  available  over  and  over  again  places 
restrictions  on  the  number  of  lattice  points  used.   With  today's 
bulk  storage  computers  this  problem  is  reduced  to  a  minimum. 

One  other  preliminary  must  be  disposed  of  before  a  problem 
of  this  type  can  be  read  into  a  computer  for  calculation.   We 
must  have  a  starting  point.   An  initial  guess  must  be  made  for 
the  values  of  <j>  at  the  interior  points.   A  typical  boundary- 
value  problem  for  elliptic  equations  is  stated  in  the  definition 
of  the  Dirichlet  problem  on  page  3.   Using  this  definition,  we 
will  prove  the  following  theorem: 

THE  MINIMUM-MAXIMUM  THEOREM:   Consider  a  harmonic  function  <t>(x,y), 
continuous  in  some  closed  bounded  region  G  =  G  +  B.   Then  the 


10 


values  of  ij)  in  G  cannot  exceed  its  maximum  on  the  boundary  B  nor 
can  they  be  less  than  its  minimum  on  B. 

PROOF:   Let  m  denote  the  maximum  of  $  on  B.   Assume  that,  at 
some  points  of  G ,  <$>  assumes  values  greater  than  m.   Then  the 
maximum  M  of  <J>  in  G  must  also  be  greater  than  m,  and  this  maxi- 
mum must  be  assumed  at  some  interior  point  Q  of  G.   We  now 
translate  the  origin  to  the  point  Q.   Under  this  transformation 
(J>  remains  harmonic.   Now  consider  the  function 

v(x,y)  i  <J>(x,y)  +  ,— (x   +  y  ), 

2d 

where  d  is  equal  to  the  least  upper  bound  of  the  distances 

between  pairs  of  points  of  B  (the  maximum  distance  across  the 

2     2     2 
region).   If  (x,y)  is  in  G,  then  x   +  y   <  d  .   From  the  above 

equation,  we  can  see  that  v(0,0)  =  $(0,0)  which  is  equal  to  M. 

On  the  other  hand,  if  a  point  (x,y)  belongs  to  the  boundary  B  of 

G,  then 

v(x,y)  <  m  +  (1/2)(M  -  m)  =  (l/2)(m  +  M)  <  M. 

This  can  be  obtained  by  noting  that  <J>(x,y)  <^  m  and  since 

2     2     2   . 
x   +  y   <  d  ,  it  follows  that 

2  4.   2 

x    v    <  1. 

Consequently,  v(x,y),  like  <(>(x,y),  must  attain  its  maximum  at  an 
interior  point  of  G.   However,  for  all  points  of  G,  we  have 
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3    v         3    v    _    3    <J>    .     3    <t>    .     2(M    -    m)         n 

2  2"  2  2  2 

3x  3y  3x  3y  d 

This  contradicts  the  fact  that  for  a  maximum,  none  of  the  second 
derivatives  of  a  function  can  be  positive . 

To  prove  the  minimum  part  of  the  theorem,  apply  the  above 
result  to  -<J>(x,y).   This  means  that  the  maximum  value  of  -<j>(x,y) 
is  the  minimum  value  of  $(x,y)  and  hence,  the  minimum  cannot 
occur  at  an  interior  point. 

The  initial  values  assigned  to  <£  could  be  zero.   However, 
the  theorem  indicates  that  the  guesses  should  lie  somewhere 
between  the  minimum  and  maximum  observed  on  the  boundary.   When 
we  intend  to  find  a  solution  for  such  a  problem  using  a  computer, 
the  accuracy  of  the  guess  is  not  of  major  significance.   At 
worst,  a  bad  guess  would  only  require  more  iterations.   After 
approximating  the  potential  at  all  interior  points  once  (one 
complete  iteration)-,  we  have  a  set  of  values  as  shown  in  table  1, 
appendix  A.   The  final  solution  for  this  example,  along  with  a 
FORTRAN  program,  is  also  illustrated  in  appendix  A. 

The  approximation  to  the  Laplacian,  as  defined  above,  is 
derived  by  neglecting  fourth  and  higher  differences  in  a 
difference  equation.    If  this  same  difference  equation  is 
considered,  neglecting  only  sixth  and  higher  differences,  we  have 
another  approximation  to  the  Laplacian.   Thus,  the  Laplacian  may 


6.   This  discussion  is  illustrated  in  detail  in  Numerical 
Analysis  by  K.  S.  Kunz,  pp.  279-80. 


be  expressed  as 


'h^OO  =  ^(6x*00  "  T7  6x*00  +  6J*00  "  17  ^OO5 
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12h 


2  [-(<i)20  +  *02  +  *-20  +  *0-2>  +  16(*10  + 


'01 


-10  +  *0-l}  -  50  *nn^  . 


>*00J 


or,  symbolically, 


h-00   ^T 


16 

16     -60     16 
16 


This  approximation  should  yield  better  results  than  the 
Laplacian  derived  for  four  neighboring  points.   It  is  very 
seldom  used,  however,  due  to  the  added  complications.   We  now 
must  know  the  potential  at  points  adjacent  to,  but  outside  the 
boundary  or  we  must  assign  weights  other  than  those  illustrated 
in  the  equation  above.   In  figure  1,  for  example,  we  do  not  have 
two  boundary  points  defined  to  the  left  of  $, . 
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LIEBMANN  METHOD  -  THREE  DIMENSIONAL 


It  would  be  suitable  to  consider  next  the  potential  when 
the  system  is  extended  to  three  coordinates.   It  turns  out  that 
for  a  simple  three-dimensional  region  it  is  not  always  possible 
to  solve  the  Dirichlet  problem.   A  restriction,  as  stated  by  I. 
G.  Petrovsky,  is  that  the  function  have  a  sufficient  number  of 
continuous  derivatives  on  the  bounding  surface.   This  condition 
can  be  stated  that  the  function  at  the  bounding  surface  be 
sufficiently  smooth. 

The  Laplacian  of  <t>  expressed  in  Cartesian  coordinates  in 
three  dimensions  is  written  as 


a2*   a2* 


hT0 


3x    3y     3z 
Therefore  we  can  use  an  approximation  similar  to  that  used  for 


(x0,y0,z0) 


Fig- 
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the  two  dimensional  case.   In  a  three  dimensional  lattice  this 
Laplacian  can  be  approximated  by  replacing  the  second  derivatives 
by  their  expressions  in  terms  of  finite  differences. 

The  Liebmann  derivation  of  a  formula  for  two  dimensions  is 
discussed  in  detail  in  chapter  12  of  Numerical  Analysis  by  K.  S. 
Kunz.   This  is  known  as  the  Laplacian  for  two  dimensions  and  is 
also  mentioned  briefly  in  the  introduction  of  this  report.   The 
following  discussion  is  analogous  to  that  of  Kunz.   The  only 
significant  difference  is  that  we  now  have  a  three  dimensional 
array. 

Suppose  the  Laplacian  is  desired  at  (xQ,yQ,zQ)  in  figure  2. 
We  can  define  u,  v  and  w  by  the  following  equations: 


x  -  xQ         y  -  yQ          z  -  zQ 
u  =  — c v  =  — , w  =  — c 


Let 


*rst  "  *(xr'yS'zt)  • 

7 

Then  by  the  use  of  Stirling's  interpolation  formula. 


'^'Vo'V  s  *ooo  +  u^x*ooo  +  (1/2!)u26x*ooo  + 


(l/3!)u(u2  -  1)U«1#000  ♦  UA!)u2(u2  -  1)«>000  ♦ 


and  hence, 


7.   Kaiser  S.  Kunz,  Numerical  Analysis ,  McGraw-Hill  Book 
Company,  Inc.,  New,  London,  Toronto,  1957,  pg.  71. 
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,2 


3  ij>(x,y0,20) 


3x 


,2 


1 


<0(x,yo,zo) 


2  2 

x  =  x_    h         3u     |  u  =  0 


"  77  "x^OOO  -  T7  6x*000  +  —) 
h 

where  the  subscript  x  on  the  <5   indicates  that  the  differences 
are  formed  with  respect  to  x. 


and 


6x*ooo  "  *-100  "^OOO  +  ""100 


6x*000  =  *-200  "  4<1>-100  +  6<|)000  -  4*100  +  *200 


Since  the  second  derivatives  with  respect  to  y  and  z  can  be 
handled  similarly,  the  following  equation  can  be  obtained. 


/s2*        32<j>        3  2<J>\  .    1    A.  2.  1      .1,  ,  \. 

(—7  +   77  +  77?  "   77l(6x*000   "   17  Sx*000   +   ~" >)  + 

\3x  3y  3z    /  x0,y0,zQ        h    \  / 

[5] 

M(6y*000  "  17  6J*000  + 1  +   ^■("z^OOO  -  T7  5z*000  +  "~ >) 

It  is  customary  to  use  only  the  approximation  to  the 
Laplacian  obtained  by  neglecting  the  fourth  and  higher  differences 
in  Eq.  5.   The  results  in  other  cases  would  be  similar.   If  the 
lattice  is  set  up  so  that  all  net  spacings  are  the  same,  h  =  j  =  k. 
Therefore 


32*    3%    32<f>\  2   2 

7T   7~2"   7T  _  Vooo 

3x     oy     3z  /  x0,yQ,zQ 
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"here  _2  _    1       ,  .2,  j.    r2  a    ,c2a         ~i 

Vh*000    "   75   (6x*000    +    6y*000    +    <5z<i)000) 

h 


72  (*100  +  ^OlO  +  *001  +  *-100  +  *0-10  +  *00-l  "  ^ooo*' 
n 


Then,  in  symbolic  notation, 


vhv000 


Using  single  subscripts  to  designate  the  relative  position  of  the 
i> '  s  in  the  net,  the  expression  can  be  written  as 


0 


vh*( 


If  one  should  choose  to  neglect  only  the  sixth  and  higher 
differences , 


12h 


-1    X^A 


6„*nnn> 


-1   *4A 


i2h2   yv°00'      i2h2   z¥°00 


would  be  evaluated. 
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Hence , 


52,    -2,         ,2. 

2      2      2 

3x    3y     oz        x0>y0,zQ 


=  V'h*000 


72  C*-100  "  2*000  +  *1C0  '  If  (<1)-200  "  ^-100 


5*000  "  4<i,ioo  +  *200)  +  *0-10  "  2*000  +  *010 


12    C*0-20 


^o-io  +  5*ooo  -  4*oio  +  *o2o)  +  *oo-a 


2*00O    +    *001    -    12    (*00-2    -    4*00-l    +    6*000    "    4<l>001 


*002)] 


12h 


2    [-C(i,-200    +    *0-20    +    *00-2    +    *200    +    *020    +    <f>002) 


+  "(*_10Q  +  *0_10  ♦  *QQ;X  ♦  *100  ♦  *010  +  *001) 


90*000] 


In  symbolic  notation, 


As  mentioned  previously,  this  approximation  is  very  seldom  used 
because  of  the  additional  values  required  at  points  adjacent  to 


but  outside  the  boundary. 

Now  consider  a  simple  problem  similar  to  the  one  discussed 
for  the  two-dimensional  case.   Assume  that  the  lattice  for  this 
example  forms  a  cube  extending  a  distance  of  4h  in  each  dimension. 
Hence,  we  have  a  lattice  of  125  points  including  27  interior 
points.   These  interior  points  can  be  numbered  in  turn  beginning 
at  the  upper  left  hand  corner  of  the  first  interior  plane.   Move 
to  the  right  to  the  last  interior  point  in  this  row.   Then  con- 
tinue to  the  second  row,  the  third  and  so  on.   When  this  plane 
is  completed,  continue  to  the  second  interior  plane. 

The  lattice  points  on  the  exterior  are  also  numbered 
beginning  at  the  upper  left.   Number  the  first  plane  by  rows 
moving  from  left  to  right.   The  planes'  containing  interior  points 
will  be  numbered  clockwise,  beginning  in  the  upper  left  hand 
corner.   The  last  plane  will  be  numbered  using  the  same  pattern 
used  on  the  first  plane.   See  figure  4  for  an  example  of  the 
numbering  of  lattice  points. 


Fie.  3 
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In  the  future,  lattice  points  will  be  referred  to  using 
single  subscript  notation.   The  single  subscript  will  refer  to 
the  point  designated  by  this  subscript  in  figure  4. 
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Assume  that  we  have  a  function  defined  on  the  boundary  which 
gives  a  maximum  potential  of  one  at  ij>32  (see  figure  3),  and  the 
values  at  other  points  of  the  boundary  as  given  in  figure  4. 
Label  the  edges  of  the  region  from  1  through  12.   Let  the  poten- 
tial at  all  points  on  edges  1,  3,  4,  5,  8,  9,  10,  11  and  12  be 
zero.   Also,  let  all  points  on  the  three  planes  determined  by 
these  edges  have  a  potential  of  zero.   The  potential  at  the 
remaining  boundary  points  can  be  defined  by  imagining  the  poten- 
tial as  the  thickness  in  the  walls  of  the  region.   Note,  at 
point  32,  the  thickness  is  one.   The  fifth  plane  is  not  shown  as 
the  potential  at  each  point  is  zero.   The  lattice  points  are 
numbered  following  the  same  pattern  used  to  number  the  first 
plane . 

To  demonstrate  the  procedure  using  the  three-dimensional 
Laplacian,  determine  the  first  approximation  to  <j>..   Since  the 
maximum  is  one  and  the  minimum  is  zero  on  the  boundary,  let  the 
initial  guess  again  by  0.5  for  all  interior  points.   From  the 
formula  derived  for  the  three-dimensional  Laplacian,  we  need  the 
values  for  six  neighboring  points.   For  $.  the  six  points  are  34, 
54  and  58  on  the  boundary  and  2,  4  and  10  from  the  interior. 
Hence , 

[5]         *x  =  1/6  <*2  +  *4  +  *1Q  ♦  034  ♦  ?S4  ♦  *B8) 

=  1/6  (.5  +  .5  +  .25  +  .5  +  .25  +  0)  =  0.333 

This  same  procedure  will  be  used  on  all  interior  points,  in  order. 
For  this  problem,  there  will  be  a  total  of  27  equations,  one  for 
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each  interior  point.   Equation  6  can  be  written  in  the  form: 

-6^  +  *2  +  *4  ♦  *1Q  =  -  (*34  ♦  *51+  ♦  *68). 

Consider  the  matrix  of  the  coefficients  of  the  system  of  equa- 
tions .   Again,  the  matrix  is  strongly  diagonal.   The  absolute 
value  of  each  diagonal  element  is  at  least  as  large  as  the  sum 
of  the  remaining  elements  in  that  row  or  column. 

When  the  values  at  all  interior  points  have  been  corrected, 
the  process  begins  over  again  at  <j.  .   When  the  value  of  the 
potential  at  all  interior  points  converges ,  to  a  specified  accur- 
acy ,  the  solution  has  been  obtained.   It  is  obvious  that  by  hand, 
this  method  would  prove  to  be  very  tedious  and  slow,  even  for 
this  relatively  "small"  problem.   However,  the  method  is  easily 
adapted  to  a  computer  and  the  FORTRAN  program  along  with  the 
results  are  given  in  appendix  B.   Note,  that  only  the  second, 
third  and  fourth  planes  of  lattice  points  are  listed,  as  the 
first  and  fifth  consist  only  of  boundary  points  with  constant 
values.   Observe  also  that  the  values  at  the  boundary  points  of 
the  planes  illustrated  remain  constant. 
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EXTENSION  OF  GAUSS-SEIDEL  METHOD 

The  Liebmann  improvement  matrix  techrr'  i    is  nothing  more 
than  the  Gauss-Seidel  method  applied  to  a  particular  problem. 
The  first  extension  of  the  Gauss-Seidel  takes  a  set  of  equations 
with  a  strongly  diagonal  matrix,  and  uses  the  latest  approxima- 
tions for  all  variables.   It  can  also  be  done  using  matrices,  but 
it  does  require  inversion  of  a  triangular  matrix. 

To  illustrate  the  method,  refer  to  the  two-dimensional 
problem  discussed  on  page  5.   From  the  first  two  difference 
equations  for  the  system,  we  can  obtain  the  following  by 
rearranging  and  solving  for  <j> . ,  where  j  is  the  index  of  the 
lattice  point: 

*x  =  1/4  ($2  +  <t>4  +  <j)11  +  *25) 

<f>2    =    1/4    (^    +    <J>3    +    4>5    +    *125     • 

Remember  that  $ . ,  j  =  10,  11,  ...  ,  25,  are  boundary  points  with 
constant  '  alues. 

Consider  the  equation  for  the  potential  4>2>  which  is  a 
linear  combination  of  <j>.  ,  <j>-,  $,.  and  <j>-,2-   Note  that  <j>,  appears 
in  the  expression  for  (j>„.   The  newest  approximation  for  <|>,  could 
be  used.   Similarly,  by  the  time  we  determine  the  potential  at 
the  third  interior  point,  we  have  a  new  value  for  cj>„.   The  value 
of  the  potential  at  all  interior  points,  except  the  first,  can 


be  expressed  in  terms  of  some  already  corrected  potential  values. 

In  this  example,  use  of  the  expression  <j>    will  represent 
the  newest  approximation  for  the  potential  of  the  internal 
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lattice  points,  s  =  1,  2,  ...  ,9.   The  approximation  immediately 
preceding  if)1    will  be  represented  by  (j>  . 

With  this  notation  in  mind,  we  can  now  rewrite  the  nine 
equations  for  this  example  as : 


>;+i 


.i+1  1 

^2  =  ¥ 

■  i+1  .  1 

ii+1  -  1 


ki+1  -  1 


.i+l 


i+1    1 

i>7    ■  4 


,i+l    1 
*8         =  ¥ 


..i+1 


*2  +  *  + 


i+1 


♦ll  +  *25) 


i  +  1 


i+1 
'2 

i+1 


+  *6  + 


si+1 


K13 


12' 


"21 


15' 
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+  <j> 


i+1 


'21 


K23 


16' 
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r 


i+1 


.i+1 
>7 

.i+1 


>9  +  *20) 


l»17  +  *19> 


Expressed  in  matrix  form,  these  equations  can  be  written  as 
6<i+1>  =  A*(i)  +  B6(i+1)  +  H 


where  K  is  a  constant  vector  derived  from  the  values  of  the 
potential  at  the  boundary  points.   The  method  is  similar  to  that 
used  to  obtain  the  matrix  coefficients  A  and  B.   These  matrices 
are  illustrated  on  the  following  two  pages. 
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This  product  is  evaluated  by  using  the  boundary  values  for 
10 ,  11 ,  . . .  ,  25 .   By  assigning  letters  to  these  matrices, 
we  can  simplify  the  appearance  of  the  expression  on  the  previous 
page. 
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(I  -  B)<j> 
(I  - 
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*<i+1)  -  t#<" 


This  matrix  T  is  an  "improvement"  matrix  for  the  Gauss-Seidel 
process.   It  is  necessary,  first,  to  calculate  T  and  G. 

After  carrying  out  the  calculations  indie.  :ed  above,  the 
improvement  matrix  T  and  the  new  constant  vector  G  appear  as 
shown  below.   Again  we  can  choose  the  initial  approximation  for 
the  interior  points  as  0.5.   For  this  example,  then 
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Compare  these  values,  using  the  Gauss-Seidel  method,  with 
the  values  shown  for  this  example  in  appendix  A,  after  the  first 
iteration.   Note,  since  the  values  using  the  improvement  matrix 
are  nearer  to  the  solution,  we  would  expect  that  this  method 
would  converge  faster  than  does  the  ordinary  Liebmann  method.   A 
FORTRAN  program  is  illustrated  in  appendix  C  and  the  results  are 
shown  for  the  1st,  5th  and  13th  iterations. 
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EXTENSION  OF  THE  GAUSS-SEIDEL  METHOD  -  THREE  DIMENSIONAL  CASE 

It  is  also  possible  to  apply  the  improvement  matrix  to  the 
three  dimensional  example  discussed  earlier.   We  will  use  the 
numbering  system  previously  described  on  page  18 .   In  the  pre- 
vious section,  we  considered  a  problem  in  two  dimensions.   It 
was  a  relatively  simple  matter  to  find  a  system  of  equations  to 
express  the  potential  at  the  interior  points.   We  then  looked  at 
the  coefficients  in  these  equations  and  placed  these  values  in 
the  proper  locations  in  one  of  the  matrices  A,  B  or  H. 

In  three  dimensions,  the  method  is  very  similar;  we  usually 
have  more  interior  points  and  larger  matrices  to  work  with.   The 
major  requirement  is  that  we  have  a  specified  pattern  for  number- 
ing the  lattice  points.   We  will  again  have  a  system  of  equations, 
one  equation  for  each  interior  point.   By  analyzing  the  coeffi- 
cients in  these  equations,  we  can  again  form  the  matrices  A,  B 
or  H  by  placing  the  values  of  these  coefficients  in  the  proper 

locations . 

st 
The  values  at  the  newest  iteration,  the  (i  +  1)   ,  can  be 

expressed  in  terms  of  the  values  at  the  i    iteration,  and  the 

matrices  A,  B  and  H  can  be  defined  as  follows: 

[7]  I«(i+1)  =  »(i+1)  =  A*(i)  ♦  B*(i+1)  ♦  H. 

s        s         s       s 

As  an  example,  consider  the  equation  for  the  potential  <j)^,  which 
is  <j>4  =  1/6  ($1    +    <t>b    +  <j>7  +  (f>13  +  <|>3g  +  <J>67). 

The  number  of  the  lattice  point  considered  determines  the  row  of 
the  matrix  in  which  the  coefficients  are  placed;  i.e.,  for  the 
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equation  above,  the  4  indicates  which  row  the  coefficients  are 
placed  in.   For  this  equation,  we  are  placing  coefficients  in 
the  fourth  row  of  the  matrices  A,  B  and  H.   Since  there  are  six 
neighboring  points  for  this  equation,  we  will  place  six  coeffi- 
cients in  the  fourth  rows  of  these  matrices.   The  matrix  in 
which  each  of  these  six  coefficients  is  placed  is  determined  by 
the  nature  of  each  of  the  six  points.   If  they  are  boundary 
points,  the  coefficients  are  placed  in  the  fourth  row  of  H.   If 
they  are  interior  points,  they  have  either  been  improved  in  this 
particular  iteration  or  they  have  not.   If  they  have  been  im- 
proved in  this  iteration,  the  index  is  less  than  the  index  of 
the  point  being  improved;  i.e.,  <t>1,    $2,  and  4>3  have  already  been 
improved  when  we  reach  <jv.   Thus,  these  coefficients  will  be 
placed  in  the  fourth  row  of  B.   If  they  have  not  been  improved 
in  this  iteration,  their  indices  are  greater  than  four.   These 
coefficients,  then,  will  be  placed  in  the  fourth  row  of  A. 

The  columns  in  which  these  six  coefficients  are  placed  are 
determined  by  the  indices  of  the  six  neighboring  points.   For 
example,  the  13  of  <J>13  implies  that  the  coefficient  of  <t>13  in 
the  expression  for  <k  will  be  placed  in  column  13.   We  already 
know  that  this  coefficient  is  also  in  row  4  of  matrix  A  by  the 
above  discussion.   Note  that  the  coefficients  of  the  six  neigh- 
boring points  are  one-sixth,  for  all  equations. 

Similarly,  this  method  car.  be  extended  to  four  or  more 
dimensions.   The  most  important  precaution  is  that  the  lattice 
Doints  be  numbered  in  order.   Secondly,  the  value  at  each  point 
must  be  improved,  in  order,  before  returning  to  the  starting 
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point.   When  the  values  at  all  interior  points  have  been  evaluated 

in  terms  of  the  six  surrounding  neighbors,  equation  7  can  be 

simplified  to 

0(i+l)  =  T^(i)  +  Q 

where  T  =  (I  -  B)_1A  and  G  =  (I  -  B)-1H. 

Substituting  the  values  for  the  boundary  points  into  the 
expression  above  and  choosing  0 . 5  as  the  original  guess  for  the 
values  at  the  interior  points,  <J>   (i  =  0)  can  be  written  as: 


(i+1) 
s 

=  0.3334 

=  0.3473 

=  0.2663 

0.4723 

0.4700 

0.3311 

0.6205 

0.5152 

0.3078 

0.3473 

0.3658 

0.2721 

0.4700 

0.4677 

0.3452 

0.5152 

0.4998 

0.3172 

0.2663 

0.2721 

0.1741 

0.3311 

0.3452 

0.2275 

0.3078 

0.3172 

0.1853 

s 

=  1,  ••-,  9; 

10,  ...,  18; 
FIRST  ITERATION 

19,  ..., 

27 

These  are  the  values  obtained  from  the  FORTRAN  program  in  appen- 
dix D.  The  matrices  A,  B,  (I  -  B)~  ,  T  and  G  are  not  shown  since 
all  but  G  have  dimensions  27  x  27,  and  would  take  up  unnecessary 
space.  If,  however,  these  matrices  should  be  required,  a  state- 
ment in  the  proper  place  in  the  FORTRAN  program  would  cause  this 
information  to  be  printed  or  punched. 


31 


EXTENSION  OF  THE  GAUSS-SEIDEL  METHOD  -  WITHOUT  MATRICES 

The  two  methods  used  up  to  this  point  are  the  Liebmarm 
method  and  an  extension  of  the  Gauss-Seidel  method.   The  Liebmann 
method  averages  the  four  neighboring  values ,  without  talcing  into 
consideration  the  newest  values.   The  extension  of  the  Gauss- 
Seidel  method  makes  use  of  the  system  of  equations  for  the  value 
of  the  potential  at  the  interior  points .   From  these  equations , 
we  manipulate  matrices  and  finally  come  up  with  an  improvement 
matrix.   Clearly,  from  the  preceding  work,  the  extension  of  the 
Gauss-Seidel  method  should  converge  much  faster  than  does  the 
Liebmann  method . 

Consider,  again,  the  two-dimensional  Liebmann  method  dis- 
cussed earlier.   Use  the  same  numbering  system  and  the  same 
pattern  for  improving  the  potential.   In  the  FORTRAN  program  for 
the  Liebmann  method  for  this  example  (appendix  A) ,  there  is  a 
statement  analogous  to 

,(k+l)  _  1/k  ,.(k)    .  ,(k)    .  .(JO    .  ,(k)   , 

9       =  1/4  L9-  -,  .  +  9.  .  .  +  d> .  .,,  +  0 .  .-.     .J 

s  l-l.:     1.3-1     1.3+1     i+l=j 

where  s  =  1 ,  2 ,  . . . ,  9 ;  i  =  row  index;  j  =  column  index.   This 

says  simply  that  at  each  interior  point  the  value  for  the  poten- 

st 
tial  of  the  (1  +  1)    iteration  is  equal  to  the  average  of  the 

values  of  the  four  neighboring  points  at  the  i    iteration. 

Notice,  that  by  altering  this  expression  only  slightly,  we 

have  another  method  for  approximating  the  value  of  the  potential. 

This  method  is  to  use  the  newest  values  as  they  are  available; 

i.e.,  by  the  time  we  correct  the  value  of  the  potential  9  ? ,  we 
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already  have  new  values  for  the  potential  at  points  1  through  6. 
Hence,  the  expression  would  now  appear  as 

*Ck+l)  _  1/1+  [(f(k+l)  +  ^Ck+1)  +  ^Ck)   +  ^(k) 
£  i-l, J    i5j-l    i>3+l    i+l, 3 

As  the  procedure  implies,  this  is  the  method  of  successive 

displacements,  sometimes  called  the  Gauss-Seidel  method.   A 

<^t    th        th 
FORTRAN  program  with  the  values  after  the  1,5'  and  13 

iterations  appears  in  appendix  E.   Compare  these  values  with 

those  shown  in  appendix  C  for  the  method  using  the  improvement 

matrix.   Note,  the  respective  values  are  exactly  the  same  in 

like  numbered  iterations. 

The  significance  of  this  observation  is  if  matrix  techniques 

are  utilized,  then  the  Gauss-Seidel  method  should  be  used  for 

computer  work.   This  method  does  not  depend  on  the  inversion  of 

matrices.   However,  if  matrices  are  not  used,  it  is  better  to 

use  the  extension  of  the  Gauss-Seidel  because  it  requires  less 

storage  then  the  Gauss-Seidel  and  converges  faster  than  does  the 

Liebmann  method. 
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UNEQUAL  ARM  METHOD 

The  methods  described  in  the  previous  pages  are  all  well 
adapted  to  calculating  the  value  of  the  potential  at  the  interior 
points  of  lattices  in  which  all  points  lie  exactly  on  the  nodes 
of  the  lattice.   However,  the  boundary  points  of  the  region  under 
observation  will,  in  general,  fall  between  lattice  points. 

This  presents  a  new  problem.   We  may  have  several  points 
that  have  one  or  more  of  their  neighboring  lattice  points  lying 
outside  the  boundary  (i.e.,  1,  3,  5,  7  and  8  in  figure  5).   The 
potential  at  the  interior  points  must  be  expressed  using  boundary 
points.   For  example,  <(>.  would  be  expressed  in  terms  of  $ .  ,  (J)^, 
<j>2  and  $4. 


1 

9"^ 

2 

a 
3\b 

4 

5 

6          1 

11 1 

7    f 

s      / 

A 

Jg 

e\ 

10 


s.h 
—} 


skh 
*1     S?h 


Fig.  5 


s4h 


f4 
Fig.  6 


Figure  6  illustrates  the  position  of  <j>,  and  its  neighbors.   The 
conceot  which  must  be  developed  here  relies  on  the  different  dis- 


tances from  <f>,  to  each  of  its  neighboring  points. 
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Before  a  solution  of  this  type  of  problem  can  be  attempted, 
a  method  must  be  developed  to  express  the  potential  at  points 

for  which  the  "arms"  are  unequal.   To  do  this  we  must  consider 

8 
the  first  and  second  divided  differences . 

f(Xl)    -   f(xQ)  fUx)  f(xQ) 

1        »■»•       *<*0»*1>    -  Xl-xQ  =x-^l^-x7^- 

f(x    ,xn)    -   f(x    ,x    ) 
2nd   D.D.       f(VXl,2)    -  x0    -    x2    ' 

f(xa)  f(xQ)  f(x2) 

"    -   <xx    -   x0)lx2    -   x^    -    (xQ    -    x2)(x1    -    xQJ    -    (x0    -    x2Ux2    -    x±) 

We   want   a   solution   for    Laplace's    equation   at    all   points    in   the 
region.      Hence,   we  must  have 

2  2 

L_i  +  ±A 

3x2         3y" 


11*    +    i'*    -    -2 


2    •    ~T-    V    *    =    0. 


If  we  approximate  our  function  <j)(x,0)  by  Newton's  divided 

9 
difference  interpolating  formula,   we  can  write 

<j>(x,0)  ■  f(x,0)  =  f(xQ)  +  (x  -  xQ)fCx0,x1)  + 


(x  -  xQ)(x  -  x1)f (xQ,x1,x2) 


Hence , 


,2 
i 

3x 


Z£^X«,X-,  ,x„j, 


8.  Kaiser  S.  Kunz ,  Numerical  Analysis ,  McGraw-Hill  Book 
Company,  Inc.,  New  York,  London,  Toronto,  1957,  chapter  5. 

9.  Ibid. 
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In  the  definition  for  f  (xQ  ,xx  ,x2 )  ,  let  x±    -    4>-,    *2    "    *2  and 
xn  =  $.  .   See  figure  6  for  an  illustration  of  the  points  $., 
i  =  1,  2,  t,  j,  k.   Note, 


*i  "  xo  =  *j  "  *i  =  "  sjh 

x1    =    ((>„  -  if.    -    s„h  -  (-s.h)  =  s„h  +  s.h 


«0  -  *2  *  *1  "  *2  "  "  S2h 


Then 

3x 


(x,o) .  :r      -^  h  h     \ 

2  \   -s.hCs  h  +  s.h)    -s2h(-s.h)    ^s^hTs^h  +  sThTJ 

2  r       2  <,<  2  *x     __lil_l 


Similarly, 


2 


3  <KQ,y)  =  1/h2 


r    2  h      2  h  +     2  ft  > 

\VSk    +    V  S4Sk         SktSk    +    VJ 


3y 

Hence , 

7  7  (  2    *4  2    *i  2    *k 

^7   +    s.(s2    +    s.)    +    sk(sk    +    s^) 


r     2  n 

ts4(sk +  s 

i2     _  p i_vA 

^^  lSjS2         sks4/   XJ 


i2I¥2 

For  an  example,  consider  a  region  defined  by  the  function 

2       2 
f(x,y)  =x/9+y/4-l=0.   This  is  a  very  specialized  sample 

problem;  however,  it  makes  use  of  the  formulas  developed  above. 

In  this  case  only  two  arms  are  shortened  for  <j>..   If  we  let  s„  = 

s„  =  1 ;  s .  =  s ;  s,  =  t  and  consider  figure  7,  equation  8  reduces 
4        j        k  ° 
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FiE.  7 


to 
[9] 


1/h' 


?K   + 


|r-T-t*6  +  s(l  +  s)*j  +  t(l  I    tT*k 

(s  +  t)A  \  _ 
~it *lj  " 


1  +  sr2 


Since  <j>  •  and  <j>,  are  values  at  boundary  points,  they  have  constant 

3       k 

values  and  equation  9  can  be  further  simplified  to 


^  =  A(J)2  +  Bcj)6  +  C 


where  A,  B  and  C  are  constants  that  are  calculated  once  and  for 
all  for  that  particular  "star". 

Assume  that  a  function  has  been  defined  so  that  the  boundary 
points  take  on  the  values  <j>„   =  <j>    =  0.25;  $    =  <j>    =  0.50; 

*31  =  *17  =  °-75;  *16  =  1-°  and  a11  other  °-°- 
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Then  the  values  of  the  unequal  "arms"  can  be  computed  from 

the  equation  for  the  curve.   For  example,  at  point  29,  f(x)  =  1. 

When  we  determine  the  "x"  we  can  compute  the  distance  s .   Using 

this  method,  then,  we  have  s  =  0.5980  and  t  =  0.4906.   Equation  9 

then  becomes  <t>,  =  0.1637  <J>0  +  0.1808  <i>   +  0.2548.   Hence,  if  our 
12b 

initial  guess  for  all  interior  points  of  the  ellipse  is  0.5,  our 
first  approximation  for  *   is  *,  =  0.1687  (.5)  +  0.1808  (.5)  + 
0.2548  =  0.4296. 

From  the  symmetry  of  the  boundary  values  for  this  example, 
we  would  expect  that  the  potential  at  the  interior  points  would 
also  be  symmetrical.   That  is,  the  potential  would  be  the  same  at 
points  symmetric  to  the  y-axis .   Table  15  of  appendix  F  illus- 
trates this. 

For  this  potential,  figure  6 


Consider  the  potential 
would  become: 


14' 


'13 


9 

s,  h 
k 

K 

s.h 

3 

*14 

s.h 

<i> 

S||h 

=    s 

*23 

15 


Then  s_  =  s,  =  s. 
2     k 


=  1: 


Fig- 


s .   From  equation  8 ,  we  have 


/    2*23   +   2^13   +    2^9     +   2*15_    /  2  +   2  \  .   1 
jTMl  +  s)    1(1+1)    1(1  +  s)    1(1  +  1)    ^1-1    l-sj^j 


This  reduces  to: 
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From  the  conditions  on  the  boundary,  we  have  that  (t>23  =  0.   Then 
the  first  term  can  be  eliminated  and  we  have 

2  2s  +  2  ,    _  n 

*i3  +  IT-IT  *9  +  *is  "— —  ht  -  °- 


It  follows  that 


*14  =  2TT~2   (*13  +  Hhf  *9  +  ♦lSj 

which  is  in  the  form 

*14  =  A*13  +  B*9  +  C*15  • 

From  the  equation  x2/9  +  y2/4  -  1  =  0 ,  we  can  find  the  value  of 

y  =  f(x,y)  =  1.3856  at  point  23.   So  s  =  1.8856  -  1.0  =  0.8856. 
Hence , 

(j>        =    0.2348<j>13    +    0.2491<j>g    +    0 .23i+8cf.15     . 

The  FORTRAN  program  for  this  example  makes  use  of  the  newest 
values  as  they  are  available.   Hence,  for  $ln,    we  have  already 
found  new  values  for  <t>g  and  <j>13 .   Combining  these  values  with 
the  initial  guess  of  0.5  in  the  first  iteration,  we  have 

A    =  0.2348^  +  0.2491cj>g  +  0.2348  CO. 5) 

=  0.2348(j>13  +  0.2491((>g  +  0.1174 

See  ap^-r.dix  F  for  the  FORTRAN  program  and  potential  values. 
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PURPOSE 

TG  DETERMINE  THE  POTENTIAL  FOH  ALL  INTERIOR  POINTS 
OF  A  TWO  DIMENSIONAL  LATTICE. 

PARAMETERS 

V   -   NUMBER  OF  POINTS  IN  THE  X  DIRECTION. 
N   -   NUMBER  OF  POINTS  IN  THE  Y  DIRECTION 
G   -   INITIAL  GUESS  FOR  THE  INTERIOR  POINTS 

REMARKS 

IN  THE  OIMENSION  STATEMENT.  FIRST  SUBSCRIPT  IS  A 
COUNTER.   THIS  COUNTS  THE  NUMBER  OF  ITERATIONS 
REQUIRED  TO  OBTAIN  THE  ACCURACY  DESIRED  BY 
STATEMENT  53. 

METHOD 

USES  THE  LIEBMANN  METHOD 


XPHI ( 30.5.5) 


DIMENSION  PHK30.5.5), 

1  FCRMAT  ( 12 ) 

2  FCRMAT  (F8. 4) 

3  FCRMAT  (5FS.2) 

1  I  FORMAT  (SF8.4//) 

1=0 

READ (  I  .  1  )M 

READ(  1.  I  )N 

READ! 1  ,2)G 

RFAD(1.3)((PHI(I.J.K),K=1,N),J=1.M) 

M l=M-l 

M=N-1 

00  2  1  J=2.M1 

DO  21  K=2,N1 

2  1  PHI ( I ,J.K)=G 

22  1=1+1 

23  IFU.GE.21)  GO  TO  42 
DO  31  J=l ,M 

DC  31  K=l ,N 
[F( J.EQ.l )GO  TO  51 
IF(K.E0.1)GD  TO  51 
IF(J.EO.M)GO  TO  51 
IF (K.EC.N)GO  TO  51 

PHI(I,J,K)=(PHI(I-1,J-1,K)+PHI(I-1,J,K+1)+PHI(I-1.J+1,K)+ 
1PHI ( 1-1 . J.K-1 ) )/4. 

31  CONTINUE 

32  WRITE (3. 11  ) (  (PHI  (  I . J.K) ,K=1 ,N) , J=l  ,M) 
GC  TO  52 

51  PH  I  <  I . J,K)=PHI ( 1-1 , J.K ) 
GO  TO  31 

52  DC  4  1  J=2. M I 
CC  41  K=2.Nl 
XPHI(I,J,K)=ABS(PHI(I,J,K)-PHI(I-1,J,K)) 

53  IF (XPHI (  I, J.K) .GE. .000 1  )GU  TO  22 
41  CONTINUE 


42 


42  STOP 


0.0000       0.5000  1.0000  0.5000  0.0000 

0.0000       0.3750  0.6250  0.3750  0.0000 

0.0000       0.3750  O.5C0O  0.3750  0.0000 

0.0000       0.2500  0.3750  0.2500  0.0000 

0.0000       0.0000  0.0C00  0.0000  0.0000 

TABLE  I  1  ST  ITERATION 

0.0000       0.5000  1.0000  0.5000  0.0000 

0.0000       0.2712  0.4383  '     0.2712  0.0000 

0.0000       0.1348  0.2031  0.1348  0.0000 

0.0000       0.0569  0.0812  0.0569  0.0000 

0.0000       0.0000  0.0000  0.0000  0.0000 

TABLE  2  10  TH  ITERATION 


0.0000 
O.COOO 
0.0000 
0.0000 
0.0000 


0.5000 
0.2636 
0. 1253 
0.0494 
0.0000 

TABLE  3 


1 .0000 
0.4289 
0.  1QQ0 
0.0717 
0 . 0000 


0.5000 
0.2636 
0. 1253 
0.0494 
0.0000 


0.0000 
0.0000 
0. 0000 
0.0000 
0. 0000 


20  TH  ITERATION 
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APPENDIX      B- 
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C  PURPOSE 

C  TO  DETERMINE  THE  PGTfcNTIAL  FOR  ALL  INTERIOR  POINTS 

C  OF  A  THREE  DIMENSIONAL  LATTICE. 

C 

C  PARAMETERS 

C  N   -   NUMBER  OF  POINTS  IN  THE  Y  DIRECTION 

C  M   -   NUMBER  OF  POINTS  IN  THE  X  DIRECTION 

C  MN  -   NUMUER  OF  POINTS  IN  THE  Z  DIRECTION 

C  G   -   INITIAL  GUESS  READ  IN 

C 

C  REMARKS 

C  IN  THE  DIMENSION  STATEMENT.  FIRST  SUBSCRIPT  IS  A 

C  COUNTER.   THIS  KEEPS  TRACK  OF  THE  NUMBER  OF 

C  ITERATIONS  REQUIRED  TO  OBTAIN  THE  ACCURACY 

C  DESIRED  BY  STATEMENT  53.   THE  REMAINING  SUBSCRIPTS 

C  ARE  FOR  THE  COMMON  X.Y.Z.  DIRECTION    IN  SPACE. 

C 

C  METHOD 

C  USES  THE  LIEBMANN  METHOD 

C 

DIMENSION  PHI (30. 5.5. S) .  XPH I ( 30. 5.5 .5 ) 

1  FQHMAT(I2) 

2  FCRMATCF8.4) 

3  FQRMAT(SF5.2) 
11  FCRMAT(SF12.4/) 

14  FORMAT! IH1. 14, 12HTH  ITERATION/)' 
1=0 

READC I. I )M 
READ( 1. 1 )N 
READ! 1. 1 )MN 

:ad(  i  ,2)g 
read! 1 .3) ( ( (phk i. j,k.mk),k=1.n) . j=1.m) ,mk=l .mn) 

1=0 

M1=M-1 

Ni=N-l 

MM=MN-l 

DO  21  J=2.Ml 

DO  21  K=2,Nl 

DO  21  MK=2.MN1 

21  PHI ( I . J,K.MK)=G 

22  1=1+1 

23  IF< I. GE. 21)00  TO  42 
DO  31  MK=1.MN 

DO  31  J=  1  ,M 
OQ  31  K=1.N 
IFIMK.EQ.l )GQ  TO  51 
:F( J.EO. 1 )G0  TO  51 
IF(K.EQ.1)G0  TO  51 
IF( J.E0.M)G0TO51 
IF(K.E0.N)GO  TO  51 
IF(MK.EQ.MN)GO  TO  51 

PHI (  I, J,K.MK)=(PHI { 1-1. J-l .K.NK1+PHI ( 1-1, J+l .K.MK1+PHK I- 
1-1 . J.K-l,KK)+PHI ( I-l. J.K+1 .MKJ+PHI I  1-1 , J.K.MK-1 1+PHI ( 1-1 . 
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2J.K, MK-t-1 ) )/6. 

31  CONTINUE 
WRITEC3, 14) I 

32  WRI TE<3, 1 t ) ( < (PHI ( I . J.K.MK) .K=l ,N) . J=l .M) .MK=l .  MN> 
GO  TO  52 

51  PHI ( I .J.K.MK )-PHl ( 1-1. J.K.MK) 
GO  TO  31 

52  00  41  J=2,M1 
DO  41  MK=2.MK1 
DO  41  K=2.N1 

XPHI ( I . J.K.MK )=ABS( PHI (I. J.K.MK)-PHI (1-1. J.K.MK)) 

53  IFCXPHI ( I. J.K.MK). GE.. 0001 >GO  TO  22 

41  CONTINUE 

42  STOP 
END 

0.0000       0.2500  0.5000       0.7500       0.7500 

0.0000       0.3333  0.5000       0.6250       0.7500 

0.0000       0.3750  0.5000       0.5000       0.5000 

0.0000       0.2917  0.3750       0.3333       0.2500 

0.0000       0.0000  0.0000       0.0000       0.0000 

2  ND  PLANE 

0.0000       0.2500  0.5000       0.5000       0.5000 

0.0000       0.3750  0.5000       0.5000       0.5000 

0.0000       0.4167  0.5000       0.5000       0.5000 

0.0000       0.3333  0.4167       0.3750       0.2500 

0.0000       0.0000  0.0000       0.0000       0.0000 

3  HO  PLANE 

0.0000       0.2500  0.2500       0.2500       0.2500 

0.0000       0.2917  0.3750       0.3333       0.2500 

0.0000       0.3333  0.4167       0.3750       0.2500 

0.0000       0.2500  0.3333       0.2917       0.2500 

0.0000       0.0000  0.0000       0.0000       0.0000 

4  TH  PLANE 

TABLE  4  1  ST  ITERATION 
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0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

0.0000 
0.0000 
0.0000 
0.0000 
0.0000 


0.2500 
0.Z037 
0.1620 
O. 1055 
0.0000 

0.2500 
0. 1620 
0. 1080 
0.0586 
0.0000 

0.2500 
0.105S 

0.0586 
0.0293 
0.0000 

TABLE 


0.5000 
0. J978 
0.3045 
0. 1620 
0.0000 

2  NO  PLANE 
0.5000 
0.3045 
0.2064 
0. 1080 
0.0000 

3  RD  PLANE 
0.2500 
0. 1620 
0. 1080 
0.0586 
0.0000 

4  TH  PLANE 
S 20  TH 


0.7500 
0.5740 
0.3978 
0.2037 
0.0000 

0.5000 
0.3978 
0.3045 
0. 1620 
0.0000 

0.2500 
0.2037 
0.1620 
0.1055 
0.0000 

ITERATION 


0.7500 
0.7500 
0.5000 
0.2500 
0.0000 

0.5000 
0.5000 
0.5000 
0.2500 
0.0000 

0.2500 
0.2500 
0.2500 
0.2500 
0.0000 


APPENDIX      C 
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US 


c  PURPOSE 

C  TO  DETERMINE  THE  POTENTIAL  FOR  ALL  INTERIOR  POINTS 

C  OF  A  TWO  DIMENSIONAL  LATTICE. 

C 

C  PARAMETERS 

c  M   -   NUMBER  OF  POINTS  IN  THE  X  DIRECTION 

C  N   -   NUMBER  OF  POINTS  IN  THE  Y  DIRECTION 

C  GU  -   INITIAL  GUESS 

C 

C  REMARKS 

C  SUBROUTINE  TO  INVERT  MATRIX  REOUIRED. 

C 

C  METHOD 

C  USES  THE  LIEBMANN  IMPROVEMENT  MATRIX. 

C 

DIMENSION  TN(9) .XPHI (20.9) .PPHI (20.9) 

DIMENSION  XID(9.9).BIB(10.9.9),BIBI(1,9.9).T(9.9).G(9) 

DIMENSION  PHI<25).A(9,9).BC9,9),H<9,25> ,ALAT{5.5) .HP(9) 

1  FORMAT! 12) 

2  F0RMATIF8.4) 

3  ~CRMAT(  16FS.2/9F5.2) 
=0RMAT(SF5.2) 

19  ?0RMAT(F16.4 > 
211  F0RMATI5X, 14. 12HTH  ITERATION/) 
READ! 1 . 1 )M 
READ!  1  .  I  )N 
READ( . ,2)GU 
MN=M*N 
C  QOUNDARY  POiNTS 

READ! I .3) (PHK JA), JA=1,MN) 

READ! 1 ,4) ( (ALAT( I  .  J  >  .  J=  1  .N  )  .  1=  I  .  M  ) 

f>N=N*N 

JMAX=(M-2) *(N-2) 

DO  22  I=1.JMAX 

00  22  J=1.JMAX 

At  I . J)=0. 

22  B( I. J )=0. 
MNAX=MN-JMAX 
DO  23  I=I.JMAX 
DO  23  J=1.MNAX 

23  H( I . J)=0. 
Ml=M-l 
N1=N-1 
MM=1 

DO  61  1=2. Ml 
DO  61  J=2,N1   l 

PHHHK  ,-ULATi  I-'l.  J  J+ALATI  -  .  J-l  )+ALAT(I  .  J  +  l  1+ALATC  I  +  t  .  J) 
1  )/4. 
IF<  {  I-D.EO.DGO  TO  62 
I  l  =  MM-{N-£ ) 
B(MM, II )=.25 
65  :F( ( J-l ).EQ. 1 >G0  TO  63 
I2=MM-l 
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H(MM. I  2)  =  . 25 

66  IFl I J+l ) .eQ.NlGO  TO  64 
I3=MM+ I 
A(MM.I3)=.2b 

67  IF( ( 1+1 ) .EQ.M)GO  TO  69 

A (MM.  I  4)=. 25 
GO  TO  61 

62  1  1=JMAX+J 
H(M«,  I  1 )=.2S 
GO  TO  6b 

63  I2  =  NN-I<-2 
H(KM, I2)=.25 
GO  TO  66 

64  I 3=JMAX+N+ [-1 
H(MM, I3)=.25 
GO  TO  67 

69  :4=JMAX+M+2*N-J-1 

H(MM. 14 )=.2S 
6  1  MM=MM+1 

JMl=JMAX+l 

J=l 

OC  101  I=1,JMAX 

HP( I )=0. 

DO  101  K=JM1,MN 

HP(  I  )=HP(  I  )+H(  I,K)*PHI (K) 

101  CONTINUE 

DO  111  1=1 . JMAX 

DO  III  J=I. JMAX 

IF( l.EO. J)  GO  TO  102 

XIOC I . J)=0. 

GO  TO  111 

102  XID< I. J)=l. 
ill  CONTINUE 

SUBTRACT  D  FROM  I 
DO  121  J=1.JMAX 
DO  121  K=1,JMAX 

121  aiac i.j.k)=xid< j,k)-b( j.k: 


CALL  MINV(BIB, JMAX.BIBI ) 

MULTIPLY  I-B  INVERSE  AND  A 

DO  131  1=1, JMAX 

DC  131  J=1,JMAX 

T( I . J )=0. 

DO  131  K=1,JMAX 

TCI.J)=T(I.J)+BIBI(1,I,K)*A(K.J) 
131  CONTINUE 

DO  141  1=1. JMAX 

GC I )=0. 

DO  141  K=1,JMAX 

G(I)=GII)+BIBI(1.I.K)*HP(K) 
141  CONTINUE 
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IC=l 

DO  151  1=1, JMAX 

151  PPMI ( IC. I )=CU 

152  DO  171  1=1. JMAX 
TNI  I  )=0. 

DC  171  K=1.JMAX 

TN( I )=TN( I  >+T< I.K)*PPHI ( :C.K ) 
171  CONTINUE 
tC=iC+l 
DO  1-31  1=1. JMAX 

191  PPHI ( IC, I >=TN( I)+G( I ) 

i  :c=ic-i 

toSUTE(3.21l>IIC 

WRITE(3,19)(PPHI<  IC.  I  ).  1  =  1.  JMAX) 

DO  192  1=1. JMAX 

XPHI ( I C.I )=AOS(PPHI ( I C.I )-PPHI( IC-I.I) ) 

IFCXPHI ( IC.I J.GE..OO01 )GO  TO  152 

192  CONTINUE 
..TCP 
END 

SUBROUTINE  M INV( AI , N I . AIBI ) 
C  FINDING  THE  INVERSE  OF  AN  NXN  MATRIX  USING  THE  L.F. 

C  METHOD. 

DIMCNSION  AU10.9.9),ai(9.9.9).CI(9.9.9).DI<9.9.9) 

DIMENSION  A  I  61 ( 1.9,9) .TKACEC 10) ;OETA(  10)  ,DETAIN( 10) 

DIMENSION  B:BI(1,9.9) . ADJ( 9.9.9) 

DO  406  1=1. NI 

DO  36  J=1.NI 

DO  36  K=1.NI 

IF(J.EQ.K)GO  TO  35 

Did  , J,K>=0. 

GO  TO  36 

35  Dill. J.K)=I. 

36  CONTINUE 
S=I 

TRACE( I  1=0. 
DO  106  J=l.NI 
DO  106  K=1.NI 

IF{ J-K ) 106,50. 106 
50  TRACEC I )=THACE( I >+AI ( I. J.K) 
106  CONTINUE 

TRACE ( I )=THACE( I )/S 

DO  206  J=l.NI 

DO  206  K=1.NI 

CI ( I , J,K)=0. 

CI ( I, J,K)=CI(I.J,K)+TRACE( I ) *D I ( I. J.K) 
206  CONTINUE 

DO  306  J=l.NI 

DO  306  K=1,NI 

Bid.  J.K)=0. 

:j  I  C  i  .  .  .  :.)=BI  (  I  .  J.K)* A  I  (  I  .  J.K) -CI  (  I  .  J.K) 
306  CONTINUE 
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00  356  J=l.NI 
00  356  K=1.NI 
AM  I  +  l  .  J.K )  =  0. 

00  356  KJ=1.NI 

'  AM I+l.J.K)=AMI  +  l. J.K)+AM 1,J,KJ)*BMI.KJ.K) 
356  CONTINUE 
406  CONTINUE 

1  =  1 

DC  S06  J--l.NI 
DO  506  K=1.NI 
ADJ ( I . J,K)=0. 
ADJ(  I,J,K)=ADJ(I.J.K)  +  ((-l.)**(NI-l)  )tBIINI-l,J.K) 

506    CONTINUE 

DETA(  I  )=( C-1.)**(NI-1>  >*TRACECNI  > 
DETAIN!  I  >  =  '../DETA(  I  ) 

DO    606    J=1.NI 

DC    606    K=1.NI 

AISI (  1, J.K)=0. 

AIBM1.J.K>  =  AIBM1.J.K)+DETAIN(  I  )*ADJ(I.J.K) 

606     CONTINUE 
RETURN 
END 


PHI(I)                  0.3750  0.2855  0.2635 

PHM2)                  0.5938  0.4509  0.4287 

PHI (3)                  0.3984  0.2746  0.2634 

PHK4)                  0.3438  0.1474  0-1251 

PHM5)                  0.4844  0.2100  0.1876 

PHI16)                  0.3457  0.1362  0.1250 

PHM7)                  0.2109  0.0603  0.0492 

PHM8)                  0.2988  0.0827  0.0715 

PHM9)                  0.1611  0.0547  0.0491 

TABLE     6  TABLE    7  TABLE    8 
l 

ITERATION           1     ST  5    TH  13    TH 
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c       purpose 

C  TU  DETERMINE  THE  POTENTIAL  FOR  ALL  INTERIOR  POINTS 

C  OF  A  THREE  DIMENSIONAL  LATTICE 

C  PARAMETERS 

C  M   -   NUMOER  OF  POINTS  IN  THE  X  DIRECTION 

C  N  .—   NUMOER  OF  POINTS  IN  THE  Y  DIRECTION 

C  MN  -   NUMOER  OF  POINTS  IN  THE  Z  DIRECTION 

C  GU  -   INITIAL  GUESS 

C 

C  REMARKS 

C  SUBROUTINE  REQUIRED  TO  INVERT  MATRIX 

C 

C  METHOD 

C  USES  THE  LIEBMANN  IMPROVEMENT  MATRIX. 

C 

DIMENSION  TN( 27 ) ,XPH I (20,27) ,PHI ( 125). XI 0(27,27) 
DIMENSION  A( 27,27) ,B( 27,27 ) ,H( 27, 125) ,ALAT(5»5.5) .HP (27) 
DIMENSION  a  10(27.27)  ,T(27.27) ,G(27),PPHI (20.27) 
DIMENSION  IP(27),V(27) 
COMMON  BIB, IP.V 

1  FORMAT; 12) 

2  FCRMAT(Fa.4) 

3  rORMAT( 16F5.2) 

4  F0RMAT(SF5.2> 
19  FORMAT( 1X,F8.4) 

211  F0RMAT(SX,I4.12HTH  ITERATION/) 
READ( 1 • I )M 
READ( 1 , 1 )N 
READ* 1 ,1 )MN 
READ( 1,2)GU 
MN2=M*N*MN 
C  BOUNDARY  POINTS 

READ( 1 ,3) (PHI ( JA). JA=1,MN2) 

READ! 1 .4 ) ( (( ALAT( I • J.K) . J=1,N) , 1=1 ,M) ,K=1 ,MN) 

JMAX=(M-2)*(N-2)*(MN-2) 

22  I=1.JMAX 
DO  22  J=l,JMAX 
A ( I , J )=0. 

22  B( I. J)=0. 

DO  23  I=1.JMAX 
DO  23  J=1.MN2 

23  HCI.J)=0. 
M1=M-1 
N1=N-1 
MN1=MN-1 

:■  .'-'  =  i 

DO  61  K=2,MNl 
DG  61  1=2, Ml 
-D  61  J=2,N1 

PHI  (MM)  =(  ALAT(  I.  J,K-1  )+ALAT(  I- 1  ,  J,  !<) +ALAT  (  I  ,  J-l  ,K)+ALAT(1 
1 . J+l , Ki+ALAT ( 1+1, J.K i+ALAT( I . J.K+1 ) )/6. 
IF( (K-l)-EQ. 1 )G0  TO  72 
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I  l=NM-(N-2)*(M-2) 
B(MM. I  1)  =  . 1667 

63  IF( ( 1-1 l.EO. 1 )G0  TO  73 

:a=KM-(N-2) 

B(MM, 12)=. 1667 

64  :F(U-1I.E0.1)GQ  TO  7* 
I3=MM-1 

H(CM. I3)=.1667 

66  IF( ( J+l ).EQ.N)GO  TO  75 
1  4  =  M  M  +  1 

A (KM, I  4)  =  . 1667 

:f( ( 1+1) .eq.m1go  to  76 

:l,  =  mm+cn-2  ) 

A(MM. IS)=.1667 

67  IF<  i  K+l  )  .EQ.KN1GO  TO  77 
I6=V'M+(N-2)*  !M-2) 
A(MM, I6>=.1667 

GO  TO  61 

72  !1=JMAX+N+N*C I-21+J 
H(MM.  t  1  )=.  1667 

GC  TO  63 

73  I2=JMAX+M*N+(K-2)*(2*N+2*(M-2> )+J 
H(MM,  I  2)  =  . 1667 

GO  TO  64 

74  ;3=j;',AX*HSNi;K-n*{2»Nt2«(«-2))-H2 
H(«K,  I  3)  =  . 1667 

GO  TO  65 

75  I4=JMAX+M*N+(K-2)*(2*N+2*(M-2> )+N+I-l 
H<MM, I  4)  =  . 1667 

GO  TO  66 

76  IS  =  j;-iAX  +  M*N+(K-2)*(2*N+2*(M-2>  )+2*N+M-l-J 
H(MM, 15)=. 1667 

GO  TO  67 

77  I6  =  JMAX  +  M*N-MK-1  ) * < 2*N+2* ( M-2 >  )  +  (  1-1  )*N+J-1 
H(MM, I6)=.1667 

6  1  .  M=MM+l 

JM1=JMAX+1 
DO  101  I=1.JMAX 
•  HP( I )=0. 

DO' 101  K=JM1,MN2 

HP(  I  )=HP(  I  )+H(  I.K>*PHI  (i<) 

101  CONTINUE 

DO  111  I=1.JMAX 
DO  111  J=l.j;-IAX 
[F(I.EQ.J)  GO  TO  102 
XID( I ,J)=0, 
GO  TO  111 

102  X :0I I , J)=l. 
Ill  CONTINUE 

SUBTRACT  B  FROM  I 
DC  121  J=..JMAX 
DO  121  K=1,JHAX 
121  BIBC J.K)=XIO( J.K)-8( J,K) 
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CALL  CRAM(JMAX.I) 

CALL  INVERS(JMAX) 
C 
C  MULTIPLY  I-H  INVERSE  AND  A 

00  131  I=1,JMAX 

DO  131  J=l,JMAX 

T(  I  .  J)=0. 

DQ  131  K=1.JMAX 

T(  I  .  J  )=T(  I  .  J  )<-DIb{  I  .K)*A(K.  J  ) 
131  CONTINUE 

00  Ml  I=1.JMAX 

C;  I  1=0. 

00  141  K=1.JMAX 

G(  I  )=G(  I  )+alO(  I.K)*HP(K) 
141  CONTINUE 

IC=1 

DO  151  I=1.JMAX 

151  PPHI ( IC. I )=GU 

152  DO  171  I=1,JMAX 
TNI ; )=0. 

00  171  K=l,.IMAX 

TN{ ! )=TN( I )  +  T{ l.K)*PPHI< IC.K) 
171  CONTINUE 
IC=IC+1 
DO  191  l=l,JMAX 

191  PPHI ( IC, I )=TN( I )+GC I) 
1 1 C= ic-l 

WRIT£(3.21D  I  IC 

'..RITE  (3,  19)  {PPHI  (IC.I).I=1.JMAX) 

00     192     I=1.JMAX 

XPHI i   IC, I )=ABS<PPHI ( IC, I )-PPHI ( IC-l.I > 1 

IFtXPHKIC.I  J.GE..0001  )G0  TD  152 

192  CONTINUE 
STOP 
END 

SUBROUTINE  CRAM(N.I) 
C       CROUT  REDUCTION  OF  AUGMENTED  MATRICES 

C       THIS  PROGRAM  PERFORMS  A  CROUT  REDUCTION  ON  A  MATRIX 
C       WITH  1=1.  THE  CROUT  REDUCTION  IS  PERFORMED  V.ITH  ROW 
C       :.\TtRCHANGES.     WITH  1=2,  THE  CROUT  REDUCTION  IS 
C       PERFORMED  WITHOUT  ROW  CHANGES. 

DIMENSION  ANC27.27). IP(27) .VC27) 

COMMON  AN.IP.V 

2240  FORMATCX  .5HP1 VOT  ,  13.  19H  I  S  LESS  THAN  l.E-08) 

2241  FORMATCIX  .SHP I VOT . I 3.7HIS  ZERO) 
GO  TO  (2200.2201),! 

2200  IDM'  -  . 

G0T02202 
22  0  1  I0MV=2 
C       .  ZDbCTION  OF  MATRIX 
2202  DO  2204  IDK=1.N 

V< IOK)=ABS{ AN( IDK. 1 ) ) 
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D02204  IDI=2,N 

IF(V(  IDK)-ABS< AN  I  I OK .  I DI ) )  ) 2203. 2204 , 2204 

2203  VC  tDK  )=AeiS(  AN(  I  DK  .  IDI  )  } 

2204  CONTINUE 


00  2222  IDK=l.N 


\ 


0ETR=-1 . 

IDKI--I0K-1 

00221*.  IDI=IDK.N 

DETPR=0.0 

IPC IDK-l ) 2208. 2208. 2206 
2200  002207  IDJ=l.IDKl 

22  0  7  DETPR  =  DETPR-i-AN<  IOI.IDJ)*AN(  10 J. IDK) 
2208  AN( 10  I. IDK )=AN( IDI . IDK)-DETPR 

GO  T0(2212.2225) . IDMV 
22  1  .J  OcTS=ABS< AN( IDI, I OK > )/V< IDI ) 

IF(DETS-DETR)2214,2214,2213 

2213  DETR=0£TS 

IP( IDK)=IOK-IDI 
GO  TO  2214 
2225  IP(IDK>-=0 

2214  CONTINUE 
I0K2=I0K-tP( IDK) 
DE7R=ANC IDK2. IDKJ/VC IDK2) 
;F(ABS(DETl-t)-l.E-08>2230.  2230.2232 

2230  WRITEO.2240  )IDK 

iF(ANUDK2.  IDK)  12232.2231.2232      ' 

2231  write. 3, 2241 )idk 

call  e;:it 

2232  VC 10K2)=V( IDK) 
V(  IDK  J--DETR 
002222  IDJ=1.N 
OETR=AN( IOK. IOJ) 

IF ( IOJ- IDK )22IS, 2215.2216 

2215  AN  I  IDK,  IDJ >  =  AN(IDK2. IDJ) 
GO  TO  2220 

2216  Dl£TPR=0.0 

IFC  IDK-l  ) 22 19. 22 19. 22  I  7 

2217  D02218  101=1. IDK1 

2218  DETPR=DETPR+ANUDK.IDI )*AN( IDI . IDJ) 

2219  AN( IDK. IDJ )=( AN< IDK2, I DJ J-DETPR ) /AN ( IDK, IDK) 

2220  IF{ IP( IDK ) )222l. 2222, 2222 

2221  AN( IDK2, IDJ)=DETR  J 

2222  CONTINUE 
RETURN 
END 

SUBROUTINE  INVERS  (N) 
C       AFTER  CALLING  THE  CRAM  SUBROUTINE  THIS  SUBROUTINE  WILL 
C       COMPUTE  THE  INVERSE  OF  THE  MATRIX  AN  AND  STORE  THE 
C       INVERSE  OF  AN  IN  AN. 

DIMENSION  AN (27, 27) . IP ( 27 ) . V ( 27 ) 

CCMMON  AN, IP. V 

DO  2280  :0K2=1.N 
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IDK  =  N+1-I0K2 

1CKI  =  IOK+1 

DO  2272  IDJ  =  IDK1.N 

V( IDJ )=AN( IDK, IDJ ) 
2272  AN< IDK. IDJ >=0.0 

AN ( IDK.  IDK >=1./AN(  I  OK.  IDK) 

IF(IDK-l)  2276.  2276,  2273 
227i  IDK3  =  IDK-1 

DO  2276  IDJ2  =  1.  IDK3 

IDJ  -     IDK3-H-IDJ2 

DdTPR  =  0-0 

IDJ1  =  IDJ  +  1 

IF(N-IDJ)  2275,  2275.  2279 
2270  DO  2274  IDI  =  IDJl.IOK 

2274  DETPR=DETPR«-AN(  IDK.  IDI  )*AN(  IDI  ,  IDJ) 
22  75  AN< IDK, I D J >=-DETPR/AN ( IDJ, IDJ) 

2276  DO  2278  IDJ=1.N 
QETPR  =0.0 

IF(N-IDK)  2278.  2278,  2284 
2284  00  2277  IDI  =  IDK1.N 

2277  QSTPR=OETPR+V(  IDI  )«AN(  IDI  ,  IDJ) 

2278  AN ( IDK,  IDJ )  =  AN ( IDK.  IDJ  J-OETPR 
2260  CONTINUE 

DO  2283  I0K2=1.N 
I0K  =  N+I-IDK2 
IDKP  =  IDK  -  IP( IDK) 
IF(IPIIDK))  22dl,  2283.  2283 
2281  DO  2282  IDI  =  I.N 
DETK=AN( IDI . IDK) 
AN{ IDI . IDK)=AN( IDI.  IDKP) 

2262  AN( IDI, IDKP)=DETR 

2263  CONTINUE 


RETURN 

END 

PHI ( 1 ) 

0.2037 

PHI ( 10 ) 

0. 1620 

PHI (2) 

0.3979 

PHI ( 1 1 ) 

0.3044 

PHI (3) 

0.S74I 

PHK  12) 

0.3978 

PHI (4) 

0.1620 

PHI (13) 

0. 1079 

PHI (5) 

0.3044 

PHI ( 14/ 

0.2062 

TABLE  9  — 


PHI (24) 

0.1619 

PHK25) 

0.0292 

PHI (26) 

0.0584 

PHI (27) 

0.1054 

TION 
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C  PURPOSE 

C  TO  DETERMINE  THE  POTENTIAL  FOR  ALL  INTERIOR  POINTS 

C  OF  A  TWO  DIMENSIONAL  LATTICE. 

C 

C  PARAMETERS 

C  M   -   NUMBER  OF  POINTS  IN  THE  X  DIRECTION. 

C  N   -   NUMf.lER  OF  POINTS  IN  THE  Y  DIRECTION 

C  G   -   INITIAL  GUESS  FOR  THE  INTERIOR  POINTS 

C 

C  REMARKS 

C  IN  THE  DIMENSION  STATEMENT.  FIRST  SUBSCRIPT  IS  A 

C  COUNTER.   THIS  COUNTS  THE  NUMBER  OF  ITERATIONS 

C  REQUIRED  TO  OBTAIN  THE  ACCURACY  DESIRED  BY 

C  STATEMENT  S3. 

C 

C  METHOD 

C  USES  THE  EXTENSION  OF  THE  GAUSS-SEIDEL  METHOD 

C 

DIMENSION  PHII2I  .5.5)  .  XPHK21.5.5) 

1  FORMAT  (12) 

2  FCRMAT  (F8.4) 

3  FORMAT  <5F5.2) 
II  FCRMAT  (5F8.4//) 

1  =  0 

READ(1.1)M 

READ< I. 1 )N 

READf 1.2)G 

READ( 1,3) < (PHK I,J.K),K=I.N).J=1,M) 

M1=M-1 

N1=N-1 

DO  21  J=2.Ml 

DO  21  K=2.N1 

21  PHI ( I.J,K)=G 

22  1=1+1 

23  IFCI.GE.2l)  GO  TO  42 
DO  31  J=1,M 

DO    3'      K=1,N 
IFCJ.EQ.DGO    TO    51 
IFIK.EQ.DGO    TO    51 
IF(J.EQ.M)GO    TO    51 
IF(K.EQ.N)GO    TO    51 

PHI  (  I. J,K)  =  (PHI( I. J-1,K)+PHI(I-S, J.K+D+PHK  1-1 , J  +  1,K)+PH 
1  I  (  I.J.K-1)  )/4. 

31  CONTINUE 

32  WRITE (3, 11 )((PHI( I . J , K ) . K= 1 , N) . J=l.M) 
GO  TO  52 

51  PHI ( I. J,K)=PHI ( 1-1 . J,K) 
GO  TO  31 

52  DO  41  J=2.M1 
DO  41  K=2.Nl 

XPHI  (  I,  J..,)=ABS(PHI  (  I, J.K) -PHI (1-1.  J.K)  ) 

53  IFIXPHI ( I. J.K) .GE..000I ) GO  TO  22 
41  CONTINUE 
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42  STOP 

ENQ 


0.0000       0.6000  1.0000  O.SOOO  0.0000 

0.0000       0.37S0  0.5938  0.3984  0.0000 

0.0000       0.3438  0.4844  0.3457  0.0000 

0.0000       0.2109  0.2988  0..611  0.0000 

0.0000       0.0000  0.0000  0.0000  0.0000 

TABLE  10  1  ST  ITERATION 

0.0000       O.SOOO  1.0000  0.5000  0.0000 

0.0000       0.2855  0.4509  0.2746  0.0000 

0.0000       0.1474  0.2100  0.1362  0.0000 

0.0000       0.0603  0.0027  0.0547  0.0000 

0.0000       0.0000  0.0000  0.0000  0.0000 

TABLE  11  5  TH  ITERATION 

0.0000       0.5000  1.0000  0.5000  0.0000 

0.0000       0.2635  0.4287  0.2634  0.0000 

0.0000       0.1251  0.1876  0.1250  0.0000 

0.0000       0.0492  0.0715  0.0491  0.0000 

0.0000       0.0000  0.0000  0.0000  0.0000 

TABLE  12  13  TH  ITERATION 
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C  PURPOSE 

C  TQ  DETERMINE  THE  POTENTIAL  FOR  ALL  INTERIOR  POINTS 

C  OF  A  TWO  DIMENSIONAL  LATTICE.   THE  BOUNDARY  POINTS 

C  OF  VHE  LATTICE  FORM  AN  ELLIPSE  AND  HENCE.  ALL 

C  POINTS  ARE  NOT  EQUIDISTANT  FROM  ALL  OTHER  POINTS. 

C  PARAMETERS 

C  M   -   NUMBER  OF  POINTS  IN  THE  X  DIRECTION. 

C  N   -   NUMBER  OF  POINTS  IN  THE  Y  DIRECTION. 

C  G   -   INITIAL  GUESS  FOR  THE  INTERIOR  POINTS. 

C 

C  REMARKS 

C  S4.S3.S2  AND  SI  REFER  TO  THE  LENGTH  OF  THE  ARMS 

C  ASSOCIATED  WITH  EACH  POINT.   THIS  IS  A  PROGRAM 

C  WRITTEN  ESPECIALLY  FOR  AN  ELLIPSE  AND  MAKES  USE  OF 

C  THE  SYMMETRIC  PROPERTIES  (I.E.  -  REGARDING  THE 

C  LENGTHS  OF  THE  'ARMS'). 

C 

C  METHOD 

C  USES  THE  EXTENSION  OF  THE  GAUSS-SEIDEL  AS  APPLIED 

C  TO  THE  PROBLEM  OF  UNEQUAL  DISTANCES  BETWEEN 

C  POINTS  (UNEQUAL  ARM-STAR). 

C 

DIMENSION  PHI (20.5.7).XPHI (20, 5. 7). S4 (5. 7), S3 (5. 7) 

DIMENSION  S2(S,7),S1(5.7).S(5.7) 

1  FORMAT! 12) 

2  FORM AT (F8. 4) 

3  F0RMAT(7F5.2) 

11  FORMAT( 1H1. 14, 12HTH  ITERATION/) 

12  F0RMAT(7F16.4//) 
READt 1 , 1 )M 
R£AD(  I. 1  IN 
READ! 1.2JG 

REAOf  !.  .3)  (  (PHI  (  I.  J.K)  ,K=1.N)  ,  J=l  .  M) 
C  Y=SQRT(4.-4.*X*X/9. ) 

C  X=SQRT(9.-9.*Y JY/4. ) 

IF(M.LE.3)G0  TO  101 

IF(N.LE.3)G0  TO  101 
C  ALL  LOOPS  THROUGH  STATEMENT  NUMBER  51  DETERMINE 

C  THE  VALUES  OF  S( I) , 1=1 , 2 , 3. 4. 

W=M 

HM=W/2.+l. 

IM  =  HM 

U  =  N 

HN=U/2.+l. 

IN=HN 

N1=N-1 

;<l  =  M-l 

IM1=IM-1 

INI=IN-1 

DO  2  1  J=2.IMl 

00  21  K=2.IN1 

KK=K-4 


63 


X  =  KK 

Y=SaRT(4.-4.*X*X/9.) 

IY=Y 

YY=IY 

S2( J.K)=Y-YV 

JJ=J-3 

"  =  J  J 

Y=-Y 

X  =  SCRT(9.-9-*Y*Y/'i.  ) 

IX=X 

XX=l  X 

IF(K.GE.3)G0  TO  22 

S3U.K)=X-XX 

GO  TO  23 

22  S3(J.K)=l. 

23  S4(J.K)=1. 
21  SI ( J.K)=1 . 

OG  31  K=IN. IN 
30  31  J=2.IM1 
S4(J.tO  =  l. 
S3C J.K )=1. 
S2< J.K  )=1. 
SI ( J.K)=1. 
31  CONTINUE 
;  N  2  =  I N  +  1 
DO  4 1  J=2.1M1 
DC  41  K=IN2.N1 
K2=K-2»(K-4) 
S4(  J.K)  =  S4CJ,K2) 
S3( J.K)=S1 ( J.K2) 
S2( J.K )=S2( J.K2) 
Sll J.K)=S3( J.K2) 

41  CONTINUE 

DO  42  J=IM, IM 
DO  42  K=2.N1 
S4( J.K)=1 . 
S3( J,K)=1. 
S2( J,K)=1. 

42  SI ( J.K)=l. 
IM2=IM+1 

'1  J=IM2.M1 
.  1  K=2.N1 
J1=J-IM 
J2=J-2«J. 

. , K)=S2( J2,K> 

S3( J.K)=S3C J2.K) 

S2( J.K)=S4{ J2,K) 

SI i J.K)=S1 ( J2.K) 

51  CONTINUE 

DO  61  J=2.M1 

CO  61  K=2.Nl 

61  PHI ( . , J.K)=G 

I  ■>  0  ■ 
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71  1=1+1 

DO  81  J=l.« 
DO  8  1  K=1.N 
IF< J.EQ.l )GO  TO  74 
[F(K.EQ.l)GO  TO  74 
IF<  J.EQ.MCO  TO  74 
:F{K.£C.N)GO  TO  74 

S(  J,K)=S4(J.K)*S2<  J.K)+S3(  J.K)*S1  <  J.K) 

PHU  I  . J.K)«Sl ( J.K)»S2( J.K) «S3( J.K) *PHI ( I . J-l .K)/<  ( S2( JtK) 
1+S4(J.K))*S(J.K))+S1 ( J.K)*S2( J.K)*S4( J,K)*PHI( I . J . K- 1) / ( ( 
2S1  ( J.K)+S3( J.K) )*S( J.K)  )+Sl( J.K)*S3< J.K)*S4( J,K)*PHI (I-l. 
3J+l.K)/< (S2( J.K)+S4l J.K) )*S( J.K) ) +S2 ( J. K ) *S3 ( J , K) *S4 ( J .K ) 
4:'PH1<I-1.J.K+1)/((S1(J.K)4-S3<J.K>>*S(J.K)) 
81  CONTINUE 

'..RITE  (3.  11)1 
62  WHITEI3. 12>< (PHI ( I. J.K) .K=1,N) . J=1.M) 

GO  TO  52 
74  PHI ( I . J.K)=PHI ( I-l. J.K) 

CO  TQ  81 
i>_  DO  9  1  J=2.M1 
DO  91  K=2.N1 

XPhI  (  I  .  J,;0=ABS(PHI  (  I  .J.K) -PHI  {  I-l  .J.K)  ) 
IF [XPHI ( I. J.K) .GE..0OO1 )GO  TO  71 
91  CONTINUE 
101  CONTINUE 
STOP 
'   END 


0.5000   0.7500   1.0000  0.7SOO   0.5000 

0.2500   0.4295   0.5457   0.6364  0.5943   0.44S4   0.2500 

0.0000   0.3574   0.4758   0.5280  0.5306   0.2690   0.0000 

0.0000   0.2160   0.3020   0.3325  0.3447   0.1941   0.0000 

0.0000   0.0000   0.0000  0.0000   0.0000 

TABLE  13  1  ST  ITERATION 
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0.5000   0.7500  I. 0000  0.7500   0.5000 

0.2500   0.3038   0.4700  0.5649  0.4671    0.3004    0.2500 

0.0000   0.1670   0.2725  0.3143  0.2693   0.1631   0.0000 

0.0000   0.0341   0.1309  0.1441  0.1290   0.0819   0.0000 

0.0000   0.0000  0.0000  0.0000   0.0000 

TABLE  14  9  TH  ITERATION 

0.5000   0.7500  1.0000  0.7500   0.5000 

0.2i      0.2980   0.4619  0.5574   0.4618   0.2979   0.2500 

0.0000   0.1605   0.2636  0.3060   0.2635   0.1604   0.0000 

0.0000   O.OS04   0.1257  0.1394  0.1257   0.0803   0.0000 

0.0000   0.0000  0.0000   0.0000   0.0000 

TABLE  15  18  TH  ITERATION 
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The  Liebmanri  method  for  solving  partial  differential  equa- 
tions r.akes  use  of  a  grid  of  points,  a  lattice.   The  potential 
at  any  interior  point  of  the  lattice  is  defined  to  be  one-fourth 
of  the  sue;  of  the  potentials  at  the  four  closest  neighboring 
points.   It  is  essential  that  the  potential  be  known  or  that  it 
can  be  measured  on  the  boundary  of  this  lattice.   For  a  simple 
two-dimensional  problem,  we  approximate  the  value  of  the  poten- 
tial at  some  interior  point.   Then  we  continue  to  approximate 
its  value  at  the  remaining  points ,  making  use  of  any  previously 
corrected  values. 

The  method  does  not  change  significantly  when  applied  to  a 
three-dimensional  lattice.   We  again  require  boundary  conditions 
to  be  defined.   Each  point  must  be  considered,  in  order,  as  we 
move  through  the  lattice  approximating  the  value  of  the  potential. 
The  only  significant  difference  is  that  we  now  have  six,  instead 
of  four,  neighboring  points  which  lie  at  a  distance  of  h  units. 
We  must  derive  an  expression  for  the  Laplacian  in  three  dimen- 
sions.  This  is  similar  to  the  two-dimensional  Laplacian  and  can 
also  be  extended  to  n  dimensions.   We  again  correct  the  value 
of  the  potential  at  each  grid  point  as  we  move  through  the 
lattice.   We  could  possibly  use  the  second  approximations  to  the 
Laplacian.   This  approximation  makes  use  of  the  eight  closest 
neighboring  points  in  the  lattice  for  the  two-dimensional  case. 
We  seldom  use  this  approximation,  however,  because  of  the  ad  jd 
complication  of  the  additional  points . 

Rather  than  correct  the  potential  at  each  interior  po_nt 


individually,  we  may  correct  all  points  simultaneously.   To  do 
this  we  use  matrix  operations.   This  " improvement ,:  matrix  is 
derived  from  a  matrix  containing  certain  coefficients.   The 
coefficients  are  taken  from  the  equations  obtained  when  applying 
the  Liebmann  method  to  the  value  at  each  point.   After  repeated 
applications  of  the  improvement  matrix,  the  potential  for  all 
points  will  converge  to  within  a  specified  accuracy.   These 
values  are  solutions  to  the  system  of  equations  which  approxi- 
mates the  partial  differential  equation. 

If  the  region  we  are  considering  can  not  be  defined  so  that 
all  boundary  points  fall  exactly  on  a  lattice  point,  we  need  to 
develop  a  new  method.   This  will  enable  us  to  approximate  the 
potential  at  the  interior  points  in  terms  of  the  four  closest 
neighbors ,  whether  or  not  they  are  lattice  points .   This  method 
is  known  as  the  "unequal  arm"  method.   As  a  result,  the  potential 
at  any  point  can  be  expressed  as  a  function  of  the  potential  at 
the  four  closest  neighbors  and  also  the  distance  between  each 
neighbor  and  the  point  whose  value  is  being  corrected. 


